Tema_4

1. MODELO DE REGRESIÓN LINEAL SIMPLE (RLS) CLÁSICO Y ROBUSTO

Utilizaremos los datos “engel” del paquete “quantreg”, datos utilizados en Koenker y Bassett (1982). Se trata de un conjunto de datos que consta de 235 observaciones sobre la renta y el gasto en alimentos de los hogares de clase trabajadora belgas. Se desea obtener un modelo que explique el gasto en alimentos en función de la renta familiar.

• Halle la correlación entre ambas variables.
• Ajuste el modelo de regresión lineal simple clásico y robusto.
• Evalúe el ajuste global del modelo y la bondad de ajuste.
• Realice la representación gráfica con la recta ajustada y los datos.
• Evalúe gráficamente los supuestos del modelo (diagnóstico) y la presencia de datos atípicos (outliers).
• En cada uno de los pasos interprete los resultados.

• En lugar de ajustar todos los modelos robustos que hemos visto, recomiendo utilizar la función “ltsReg” del paquete “rrcov” que permite obtener un resumen del modelo comparable al de la función “lm”.

• Comparar los resultados de los modelos lm y ltsReg.

library(quantreg)
data("engel")
head(engel)

##     income  foodexp
## 1 420.1577 255.8394
## 2 541.4117 310.9587
## 3 901.1575 485.6800
## 4 639.0802 402.9974
## 5 750.8756 495.5608
## 6 945.7989 633.7978

attach(engel)

Estudio de normalidad y graficado

boxplot(engel)

# Se aprecia que no siguen normalidad, y existen outliers en ambas variables.

shapiro.test(engel$income)

## 
##  Shapiro-Wilk normality test
## 
## data:  engel$income
## W = 0.79204, p-value < 2.2e-16

# p < 0.05. Se rechaza la hipótesis nula de existencia de normalidad.

shapiro.test(engel$foodexp)

## 
##  Shapiro-Wilk normality test
## 
## data:  engel$foodexp
## W = 0.87471, p-value = 5.45e-13

# p < 0.05. Se rechaza la hipótesis nula de existencia de normalidad.

Hallar la correlación entre variables. Se usará correlaciones robustas al existir outliers.

plot(income, foodexp, pch=20)
abline(coef(lm(foodexp~income)), col='green')

En principio sí parece haber linealidad.

library(WRS)

pbcor(income, foodexp)

## $cor
## [1] 0.9342118
## 
## $test
## [1] 39.9758
## 
## $siglevel
## [1] 0
## 
## $n
## [1] 235

wincor(income, foodexp)

## $cor
## [1] 0.9273018
## 
## $cov
## [1] 34445.98
## 
## $siglevel
## [1] 0
## 
## $n
## [1] 235

La correlación entre ambas variables es muy fuerte y positiva. p-valor menor a 0.05 por lo que rechazamos la hipótesis nula (correlación es 0). Existe una relación significativa entre ambas variables.

Ajuste del modelo regresión lineal simple clásico y robusta.

# Alimentos (VD) en función de la renta (VI). Es decir, foodexp en función de income.

summary(fitLS <- lm(foodexp~income))

## 
## Call:
## lm(formula = foodexp ~ income)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -725.70  -60.24   -4.32   53.41  515.77 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 147.47539   15.95708   9.242   <2e-16 ***
## income        0.48518    0.01437  33.772   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 114.1 on 233 degrees of freedom
## Multiple R-squared:  0.8304, Adjusted R-squared:  0.8296 
## F-statistic:  1141 on 1 and 233 DF,  p-value: < 2.2e-16

Como el intercepto (beta 0) e income(beta1), poseen un p-valor < 0.05 rechazamos la hipótesis nula por lo que ambas son significativos, deben aparecer en la ecuación.

La ecuación quedaría como: foodexp = 0.48518 * income + 147.47539

F-statistic (ANOVA), el p-value es < 0.05 nos aporta que el modelo es significativo.

Los parámetros del modelo intercepto y pendiente son ambos significativos (p-valor < 0.001)

Ajuste global del Modelo y Bondad de Ajuste…

Se nos muestra en F-statistic del paso anterior. Seria: F(1,233)=1141, p-value < 0.001

La bondad de ajuste viene dada por R2 que en este caso es de 0.8304, el modelo es muy eficaz (entre 0 y 1).

# graficamos...
plot(foodexp~income, pch=20)
abline(fitLS)

La linealidad se puede apreciar en el gráfico de los datos anterior. Existe linealidad.

Evalúe gráficamente los supuestos del modelo (diagnóstico) y la presencia de datos atípicos (outliers)

# Gráfico genérico del modelo:

par(mfrow=c(2,2))
plot(fitLS)

Independencia de residuos

library(lmtest)
dwtest(fitLS)

## 
##  Durbin-Watson test
## 
## data:  fitLS
## DW = 1.4108, p-value = 2.677e-06
## alternative hypothesis: true autocorrelation is greater than 0

Se usó prueba de Durbin-Test. Se obtuvo: Autocorrelación positiva dwtest < 2.

Vemos los gráficos anteriores uno a uno…

# Residuals vs. Fitted.
par(mfrow=c(1,1))
plot(fitLS, which=1)

# No aparece ningún tipo de patrón. Parece media 0. No hay heterocesdasticidad ni problemas de linealidad. No hay correlacionalidad en los residuos.
# Se aprecian 3 outliers 59, 105 y 138.

# Normal Q-Q
plot(fitLS, which=2)

# A excepción de los outliers comentados anteriormente, parece que se ajustan los datos a la recta.
# La hipótesis de que los errores siguen una distribución normal parece cumplirse.

# Scale - Location
plot(fitLS, which=3)

# Como en el caso del gráfico primero, no aparece ningún tipo de patrón.

# Residuals vs. Leverage
plot(fitLS, which=4)

# Muestra los 3 outliers que se deben estudiar.

# graficamos distancias de cook
cooks <- cooks.distance(fitLS)
hat <- lm.influence(fitLS)$hat

par(mfrow = c(2,1))
plot(income,foodexp, col=ifelse(cooks > quantile(cooks, .90), 'red', 'black'), pch=20)
plot(income,foodexp, col=ifelse(hat > quantile(hat, .90), 'green', 'black'), pch=20)

par(mfrow = c(2,1))
plot(income,predict(fitLS), col=ifelse(cooks > quantile(cooks, .90), 'red', 'black'), pch=20)
plot(income,predict(fitLS), col=ifelse(hat > quantile(hat, .90), 'green', 'black'), pch=20)

# Parece que hay outliers que influyen en la recta. Se debería trabajar mediante el uso de técnicas robustas (Ya comentado al inicio).

# Uso de validación cruzada para ver el ajuste del modelo.
library(caret)
datos <- engel
train(foodexp~income, data = datos, method = "lm")

## Linear Regression 
## 
## 235 samples
##   1 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 235, 235, 235, 235, 235, 235, ... 
## Resampling results:
## 
##   RMSE      Rsquared 
##   124.3693  0.8419717
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
## 

# R2 es 0.84. Posee un buen ajuste.

---

AJUSTE USANDO ltsReg PARA ROBUSTOS (Regresión por mínimos cuadrados recortados)

# Graficamos los outliers
library(car)
plot(income, foodexp, pch=20, main = "Regresión Paramétrica")
abline(coef(lm(foodexp~income)), col='green')
out <- influence.measures(fitLS)
showLabels(engel$foodexp, engel$income, row.names(engel), 
           id.method=row.names(engel)[row.names(engel) %in% names(which(apply(out$is.inf, 1, any)))])

##  49  59  61  85 105 119 125 128 138 155 158 220 
##  49  59  61  85 105 119 125 128 138 155 158 220

# Existencia de outliers. Lo conveniente es aplicar regresión robusta.

fitLSR <- ltsReg(income, foodexp)
summary(fitLSR)

## 
## Call:
## ltsReg.default(x = income, y = foodexp)
## 
## Residuals (from reweighted LS):
##     Min      1Q  Median      3Q     Max 
## -188.51  -36.66    0.00   34.26  161.44 
## 
## Coefficients:
##           Estimate Std. Error t value Pr(>|t|)    
## Intercept -5.25412   12.19031  -0.431    0.667    
## income     0.68244    0.01261  54.133   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 69.49 on 209 degrees of freedom
## Multiple R-Squared: 0.9334,  Adjusted R-squared: 0.9331 
## F-statistic:  2930 on 1 and 209 DF,  p-value: < 2.2e-16

Como el income(beta1), poseen un p-valor < 0.05 rechazamos la hipótesis nula por lo que solo él debe aparecer en la ecuación. El intercepto no pues aceptamos la hipótesis nula.

La ecuación quedaría como: foodexp = 0.68244 * income

F-statistic (ANOVA), el p-value es < 0.001 nos aporta que el modelo es significativo.

Ajuste global del Modelo y Bondad de Ajuste…

Se nos muestra en F-statistic del paso anterior. Seria: F(1,209)=2930, p-value < 0.001

La bondad de ajuste viene dada por R2 que en este caso es de 0.9334, el modelo es muy eficaz (entre 0 y 1).

# graficamos...
par(mfrow=c(1,2))
plot(income, foodexp, pch=20, main = "Regresión Paramétrica")
abline(coef(lm(foodexp~income)), col='green')
plot(foodexp~income, pch=20, main = "Regresión Robusta")
abline(fitLSR, col='red')

¿Qué curva elegir?

# Obtenemos los valores...
summary(fitLS)$r.squared

## [1] 0.8303646

summary(fitLS)$coefficients[,4]

##  (Intercept)       income 
## 1.573614e-17 9.918947e-92

summary(fitLSR)$r.squared

## [1] 0.9334254

summary(fitLSR)$coefficients[,4]

##     Intercept        income 
##  6.669065e-01 6.193203e-125

La regresión robusta obtuvo un mejor valor RSaquare, pasando del 0.83 a 0.93 con lo que ofrece un mayor ajuste.

---

---

2. MODELO DE REGRESIÓN LINEAL MÚLTIPLE (RLM) CLÁSICO Y ROBUSTO

Utilizaremos los datos “chicago”, del paquete “faraway”. Corresponden a datos de un estudio sobre la disponibilidad de seguros en Chicago desde diciembre 1977 a febrero 1978. Los datos están dados para cada código postal de Chicago (nombres de las filas o casos). Las variables son:

• race: la composición racial en porcentaje minoritario
• fire: incendios por cada 100 unidades de vivienda
• theft: robo por cada 1000 habitantes
• age: porcentaje de las viviendas construidas antes de 1939
• volact: nuevas políticas de vivienda, además de las renovaciones menos cancelaciones y no renovaciones por cada 100 unidades de vivienda
• involact: nuevas políticas del plan FAIR y renovaciones por cada 100 unidades de vivienda
• income: ingreso medio familiar

Queremos explicar la variable “involact” en función de las demás variables, excepto “volact”, y con “income” en logaritmo.

• Halle la matriz de correlaciones para todas las variables y el gráfico de dispersión. (X)
• Ajuste el modelo de regresión lineal múltiple clásico y robusto. (X)
• Evalúe el ajuste global del modelo y la bondad de ajuste.
• Estime la importancia relativa de los predictores.
• Realice la representación gráfica con la recta ajustada y los datos.
• Evalúe gráficamente los supuestos del modelo (diagnóstico) y la presencia de datos atípicos (outliers).
• En cada uno de los pasos interprete los resultados. (X)

• Nuevamente, en lugar de ajustar todos los modelos robustos que hemos visto, recomiendo utilizar la función “ltsReg” del paquete “rrcov” que permite obtener un resúmen del modelo comparable al de la función “lm”.

• Comparar entre los modelos lm y ltsReg la significación de los coeficientes, el ajuste global y la bondad de ajuste.

# Obtenemos los valores...
library(faraway)
data("chicago")
head(chicago)

##       race fire theft  age volact involact income
## 60626 10.0  6.2    29 60.4    5.3      0.0  11744
## 60640 22.2  9.5    44 76.5    3.1      0.1   9323
## 60613 19.6 10.5    36 73.5    4.8      1.2   9948
## 60657 17.3  7.7    37 66.9    5.7      0.5  10656
## 60614 24.5  8.6    53 81.4    5.9      0.7   9730
## 60610 54.0 34.1    68 52.6    4.0      0.3   8231

attach(chicago)
  
# preparación de los datos. Aplicamos logaritmos sobre la columna income, eliminamos la columna no usada en la regresión.
  
datos <- chicago
datos$income <- log(datos$income)
datos <- datos[,-5]
head(datos)

##       race fire theft  age involact   income
## 60626 10.0  6.2    29 60.4      0.0 9.371098
## 60640 22.2  9.5    44 76.5      0.1 9.140240
## 60613 19.6 10.5    36 73.5      1.2 9.205127
## 60657 17.3  7.7    37 66.9      0.5 9.273878
## 60614 24.5  8.6    53 81.4      0.7 9.182969
## 60610 54.0 34.1    68 52.6      0.3 9.015663

Halle la matriz de correlaciones para todas las variables y el gráfico de dispersión.

# graficamos a nivel general...
  
pairs(datos, pch=20)

# graficamos con correlación paramétrica
  
library(corrplot)
par(mfrow = c(1,1))
corrplot(cor(datos), type = "lower")

corrplot(cor(datos), method = "number", type = "lower")

# Gráfico resumen con distribucion (diagonal), linea de ajuste y scatter plot bivariable (inferior) y valor de correlacion y significancia (superior).
  
library(PerformanceAnalytics)
chart.Correlation(datos, histogram=TRUE, pch=19)

En el gráfico se aprecia una correlación fuerte y positiva entre involact y fire (0.70), involact y race (0.71) y poco significativa pero positiva entre involact y age (0.48)

Una correlacion fuerte y negativa entre involact e income (-0.70).

Como se aprecian outliers, realizamos la correlación robusta (en esta caso winsorizada).

library(WRS)
winall(datos)

## $cor
##            [,1]       [,2]       [,3]       [,4]       [,5]       [,6]
## [1,]  1.0000000  0.7521896  0.5300843  0.2588940  0.7305137 -0.7818053
## [2,]  0.7521896  1.0000000  0.4285737  0.4674647  0.8580696 -0.8155316
## [3,]  0.5300843  0.4285737  1.0000000  0.3386472  0.4067537 -0.4810901
## [4,]  0.2588940  0.4674647  0.3386472  1.0000000  0.5294211 -0.5667558
## [5,]  0.7305137  0.8580696  0.4067537  0.5294211  1.0000000 -0.7360743
## [6,] -0.7818053 -0.8155316 -0.4810901 -0.5667558 -0.7360743  1.0000000
## 
## $cov
##            [,1]        [,2]       [,3]       [,4]        [,5]        [,6]
## [1,] 678.375060 101.9276781 118.657817  91.087831  8.51273358 -3.25362950
## [2,] 101.927678  27.0682239  19.163367  32.853534  1.99736818 -0.67796175
## [3,] 118.657817  19.1633673  73.864015  39.315819  1.56406105 -0.66065901
## [4,]  91.087831  32.8535338  39.315819 182.476438  3.19970860 -1.22330262
## [5,]   8.512734   1.9973682   1.564061   3.199709  0.20017576 -0.05262134
## [6,]  -3.253630  -0.6779618  -0.660659  -1.223303 -0.05262134  0.02553107
## 
## $p.values
##              [,1]         [,2]         [,3]         [,4]         [,5]
## [1,]           NA 3.092819e-08 0.0002647333 0.0833663910 1.021741e-07
## [2,] 3.092819e-08           NA 0.0036608185 0.0014461527 1.162315e-11
## [3,] 2.647333e-04 3.660818e-03           NA 0.0228120339 5.934029e-03
## [4,] 8.336639e-02 1.446153e-03 0.0228120339           NA 2.699239e-04
## [5,] 1.021741e-07 1.162315e-11 0.0059340287 0.0002699239           NA
## [6,] 5.057372e-09 4.676761e-10 0.0010219429 0.0000857958 7.590220e-08
##              [,6]
## [1,] 5.057372e-09
## [2,] 4.676761e-10
## [3,] 1.021943e-03
## [4,] 8.579580e-05
## [5,] 7.590220e-08
## [6,]           NA

Se aprecian datos ligeramente superiores aunque siguen el mismo patrón. También en el caso de involact con income es fuerte y neamtiva (-0.73)

Se aprecia también mirando p-value, que en todos los casos, p-value < 0.05, por lo que rechazamos la hipótesis nula (correlación es 0). Existe una relación significativa entre las variables.

Ajuste el modelo de regresión lineal múltiple clásico y robusto.

summary(fitLM0 <- lm(involact ~ ., data = datos))

## 
## Call:
## lm(formula = involact ~ ., data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.85393 -0.16922 -0.03088  0.17890  0.81228 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.573976   3.857292  -0.927 0.359582    
## race         0.009502   0.002490   3.817 0.000449 ***
## fire         0.039856   0.008766   4.547 4.76e-05 ***
## theft       -0.010295   0.002818  -3.653 0.000728 ***
## age          0.008336   0.002744   3.038 0.004134 ** 
## income       0.345762   0.400123   0.864 0.392540    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3345 on 41 degrees of freedom
## Multiple R-squared:  0.7517, Adjusted R-squared:  0.7214 
## F-statistic: 24.83 on 5 and 41 DF,  p-value: 2.009e-11

El intercepto posee un p-value > 0.05 por lo que lo excluimos al no ser significativo. Volvemos a ajustar…

summary(fitLM <- lm(involact ~ .-1, data = datos))

## 
## Call:
## lm(formula = involact ~ . - 1, data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.87307 -0.14861 -0.01933  0.19550  0.81707 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)    
## race    0.007994   0.001881   4.250 0.000116 ***
## fire    0.036445   0.007942   4.589    4e-05 ***
## theft  -0.009563   0.002700  -3.541 0.000990 ***
## age     0.007065   0.002373   2.977 0.004810 ** 
## income -0.024709   0.015067  -1.640 0.108480    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.334 on 42 degrees of freedom
## Multiple R-squared:  0.8708, Adjusted R-squared:  0.8554 
## F-statistic:  56.6 on 5 and 42 DF,  p-value: < 2.2e-16

Ahora todas las variables son significativas, a excepción de income con p-value > 0.05. Pero le añadimos a la ecuación.

La ecuación quedaría como: involact = 0.0095 * race + 0.0398 * fire - 0.010 * theft + 0.008 * age - 1.64 * income

F-statistic (ANOVA), el p-value es < 0.001 nos aporta que el modelo es significativo.

Se nos muestra en F-statistic del paso anterior. Seria: F(2,42) = 101.2, p-value < 0.001

La bondad de ajuste viene dada por R2 que en este caso es de 0.87, el modelo es muy eficaz (entre 0 y 1). El modelo es significativo y bueno.

anova(fitLM0, fitLM)

## Analysis of Variance Table
## 
## Model 1: involact ~ race + fire + theft + age + income
## Model 2: involact ~ (race + fire + theft + age + income) - 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     41 4.5882                           
## 2     42 4.6843 -1 -0.096073 0.8585 0.3596

# p-valor > 0.05 nos dice que eliminar el intercepto no tiene efectos significativos en el modelo.
  
avPlots(fitLM)

Colinealidad:

vif(fitLM)

##     race     fire    theft      age   income 
## 3.374141 6.258661 4.712647 9.820465 8.185995

# Como los datos que arroja la prueba no son mayores de 10, nos dice que los predictores NO están correlacionados. No hay colinealidad.

Seleccion automática de predictores

summary(step(fitLM, action = F, direction = "forward"))

## Start:  AIC=-98.38
## involact ~ (race + fire + theft + age + income) - 1

## 
## Call:
## lm(formula = involact ~ (race + fire + theft + age + income) - 
##     1, data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.87307 -0.14861 -0.01933  0.19550  0.81707 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)    
## race    0.007994   0.001881   4.250 0.000116 ***
## fire    0.036445   0.007942   4.589    4e-05 ***
## theft  -0.009563   0.002700  -3.541 0.000990 ***
## age     0.007065   0.002373   2.977 0.004810 ** 
## income -0.024709   0.015067  -1.640 0.108480    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.334 on 42 degrees of freedom
## Multiple R-squared:  0.8708, Adjusted R-squared:  0.8554 
## F-statistic:  56.6 on 5 and 42 DF,  p-value: < 2.2e-16

# En esta caso las variables se incluyen en el modelo final.

Importancia relativa de los predictores.

library(relaimpo)
(crlm <- calc.relimp(fitLM0, type = c("lmg", "last", "first", "pratt"), rela=T))

## Response variable: involact 
## Total response variance: 0.4017299 
## Analysis based on 47 observations 
## 
## 5 Regressors: 
## race fire theft age income 
## Proportion of variance explained by model: 75.17%
## Metrics are normalized to sum to 100% (rela=TRUE). 
## 
## Relative importance metrics: 
## 
##               lmg       last      first       pratt
## race   0.29243355 0.24877115 0.29140971  0.46388154
## fire   0.30950866 0.35299237 0.28272657  0.54707074
## theft  0.05749694 0.22790620 0.01280705 -0.07206704
## age    0.12207995 0.15757884 0.12945724  0.18789031
## income 0.21848090 0.01275145 0.28359943 -0.12677556
## 
## Average coefficients for different model sizes: 
## 
##                  1X          2Xs          3Xs          4Xs          5Xs
## race    0.013882353  0.010877318  0.009202480  0.008628879  0.009502223
## fire    0.047902503  0.040131301  0.035927286  0.035527385  0.039856040
## theft   0.004254602 -0.002893112 -0.005717252 -0.007625214 -0.010294505
## age     0.013356717  0.008219986  0.006461752  0.006660529  0.008335600
## income -1.798818702 -1.363378514 -0.889943051 -0.344752568  0.345761521

plot(crlm)

# En la mayoría de los métodos se aprecia como "race" y "fire" los que más importancia tienen y de ellos "fire" el más imporante. (que además son las de mayor significancia con ***)

Diagnóstico

par(mfrow = c(2,2))
plot(fitLM)

Resiuals vs. Fitted.

par(mfrow=c(1,1))
plot(fitLM, which=1)

# Están bastante cercanos a 0, no hay casi dispersion, aparecen 3 outliers.

Normal Q-Q

plot(fitLM, which=2)

# Se aprecia que se ajustan lso residuos a la curva por lo que existe normalidad. Se muestran los 3 outliers

Scale-Location

plot(fitLM, which=3)

# No se aprecian enfecto estilo embudo.

Resiuals vs. Leverage

plot(fitLM, which=4)

# Se muestra la existencia de los 3 outliers.

Test de datos atípicos

outlierTest(fitLM)

## 
## No Studentized residuals with Bonferonni p < 0.05
## Largest |rstudent|:
##        rstudent unadjusted p-value Bonferonni p
## 60610 -3.262772          0.0022278       0.1047

# Solo está detectando uno, el más alejado el 60610.

Evaluación de observaciones influyentes.

avPlots(fitLM)

influencePlot(fitLM, id.method = "noteworthy", main="Grafico Influencia")

##          StudRes       Hat      CookD
## 60610 -3.2627716 0.2105705 0.46185190
## 60607  0.4509906 0.6383261 0.07318241

# Los más influyentes son los outliers 60610 y 60607 aunque con el outlierTest solo detectó el 60610.

Normalidad de los residuos

qqPlot(fitLM)

# Ya comentado anteriormente que se ajusta bastante a la recta.

Distribución de los residuos estunderizados

library(MASS)
sresid <- studres(fitLM)
hist(sresid, freq = F, main = "Distribucion de los residuos")
xfitLM <- seq(min(sresid), max(sresid), length=40)
yfitLM <- dnorm(xfitLM)
lines(xfitLM, yfitLM)

# Distribución bastante normal, aunque ligeramente sesgada a izda.

Evaluación de homocedasticidad

ncvTest(fitLM)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 9.205779    Df = 1     p = 0.002412523

# Se rechaza Ho por lo que no hay varianza constante en los errores.
spreadLevelPlot(fitLM)

## 
## Suggested power transformation:  0.4330704

Evaluación de la linealidad

crPlots(fitLM)

# Race, Fire y age son lo que más se mantienen en la linealidad.

Evaluación de la independencia de errores

durbinWatsonTest(fitLM)

##  lag Autocorrelation D-W Statistic p-value
##    1      -0.1106017      2.210188   0.552
##  Alternative hypothesis: rho != 0

# p-valor al límite 0.05.

Por existencia de outliers, se debería trabajar aplicando técnicas robustas. Que es el siguiente paso.

---

AJUSTE USANDO ltsReg PARA ROBUSTOS (Regresión por mínimos cuadrados recortados)

Graficamos los outliers

pairs(datos)

out <- influence.measures(fitLM)
plot(fitLM, which=4)

# Existencia de outliers. Lo conveniente es aplicar regresión robusta.

# library(rrcov)
fitLMR0 <- ltsReg(involact ~ race + fire + theft + age + income)
summary(fitLMR0)

## 
## Call:
## ltsReg.formula(formula = involact ~ race + fire + theft + age + 
##     income)
## 
## Residuals (from reweighted LS):
##       Min        1Q    Median        3Q       Max 
## -0.387900 -0.077283 -0.005965  0.115227  0.496925 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## Intercept -4.324e-01  3.313e-01  -1.305   0.2002    
## race       5.678e-03  1.694e-03   3.351   0.0019 ** 
## fire       4.786e-02  5.861e-03   8.166 1.03e-09 ***
## theft     -8.572e-03  1.919e-03  -4.467 7.57e-05 ***
## age        4.202e-03  1.914e-03   2.195   0.0347 *  
## income     2.025e-05  2.115e-05   0.958   0.3445    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2202 on 36 degrees of freedom
## Multiple R-Squared: 0.8475,  Adjusted R-squared: 0.8263 
## F-statistic:    40 on 5 and 36 DF,  p-value: 1.005e-13

# Todos menos income, poseen un p-valor < 0.05 rechazamos la hipótesis nula por lo que debeb aparecer en la ecuación. El intercepto no pues aceptamos la hipótesis nula.

fitLMR <- ltsReg(involact ~ race + fire + theft + age + income - 1)
summary(fitLMR)

## 
## Call:
## ltsReg.formula(formula = involact ~ race + fire + theft + age + 
##     income - 1)
## 
## Residuals (from reweighted LS):
##      Min       1Q   Median       3Q      Max 
## -0.39796 -0.04698  0.00000  0.10471  0.52330 
## 
## Coefficients:
##          Estimate Std. Error t value Pr(>|t|)    
## race    4.580e-03  1.417e-03   3.233 0.002577 ** 
## fire    4.530e-02  5.858e-03   7.732 3.07e-09 ***
## theft  -7.510e-03  1.899e-03  -3.955 0.000333 ***
## age     1.983e-03  1.463e-03   1.356 0.183368    
## income -5.939e-06  6.326e-06  -0.939 0.353900    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2263 on 37 degrees of freedom
## Multiple R-Squared: 0.9072,  Adjusted R-squared: 0.8947 
## F-statistic: 72.38 on 5 and 37 DF,  p-value: < 2.2e-16

Ahora todas las variables son significativas, a excepción de income con p-value > 0.05. Pero le añadimos a la ecuación.

Respecto al modelo paramétrico comentar que da menor significancia a la variable “race” e income no aparece con su estimado en netativo.

La ecuación quedaría como: involact = 0.00458 * race + 0.0453 * fire - 0.00751 * theft + 0.0019 * age - 0.0000059 * income

F-statistic (ANOVA), el p-value es < 0.001 nos aporta que el modelo es significativo.

Ajuste global del Modelo y Bondad de Ajuste…

Se nos muestra en F-statistic del paso anterior. Seria: F(5,37)=72.38, p-value < 0.001

La bondad de ajuste viene dada por R2 que en este caso es de 0.90, el modelo es muy eficaz (entre 0 y 1). El modelo es significativo y bueno.

Comentar que R2 el ligeramente inferior al obtenido por el anterior (paramétrico)

Estimacion de la precisión de los modelos.

La regresión robusta obtuvo un mejor valor RSaquare, pasando del 0.87 a 0.90 con lo que ofrece un mayor ajuste.

detach(chicago)
rm(list = ls())